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Abstract: We consider several stochastic service systems, and study the asymptotic behavior of the moments 

of various quantities that have application to models for random interval graphs and algorithms for searching 
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1. Introduction 

Our goal in this paper is the study of stochastic service systems, with an eye to two applications: models 
for random interval graphs and the analysis of search algorithms. The systems we study are traditionally 
designated M/M/l (independent exponentially distributed interarrival times, independent exponentially dis- 
tributed service times and a single server) and M/M/oo (independent exponentially distributed interarrival 
times, independent exponentially distributed service times and infinitely many servers). 

In a previous paper, Pippengcr [P] introduced models for random interval graphs based on stochastic 
service systems, and analyzed, among others, a model based on the M/M/oo system. For this system, 
idle periods (intervals of time during which all servers are idle) alternate with busy periods (intervals of 
time during which at least one server is busy). A random graph can be constructed by considering a busy 
period, letting the vertices correspond to customers served during this busy period, and putting an edge 
between two vertices if the service intervals of the corresponding customers overlap. Since edges correspond 
to intersections between intervals in a totally ordered set, the resulting graph is an interval graph, so this 
procedure yields a model for random interval graphs. (Other models for random interval graphs have been 
considered by Scheinerman [SI, S2] and by Godehardt and Jaworski [G2].) 

Suppose that in the system M/M/oo customers arrive at rate A and are served at rate 1. Pippenger 
[P] showed that the number N of vertices in the corresponding random graph (which corresponds to the 
number of customers served during the busy period) is such that the distribution of N/e x tends to that of 
an exponential with mean 1 as A — > oo. Furthermore, the chromatic number K of the graph (which for an 
interval graph equals the number of vertices in the largest clique of the graph, and corresponds to the largest 
number of customers simultaneously in the system during the busy period) is such that K/eX tends to 1 in 
probability as A — > oo. 

Our first goal in this paper is to study the corresponding random interval graph model for the M/M/l 
system. In this case we must have A < 1 to ensure that the busy period is finite with probability one, and 
that the number of customers served during the busy period has finite expectation, and we shall be interested 
in asymptotics as A — > 1. When there is only one server, customers who arrive when the server is busy must 
wait for service, and a service discipline (which determines which of the waiting customers will be served 
next when the current service interval ends) must be specified. As a result, the corresponding interval graph 
will depend on the service discipline used. Consider, for example, a busy period consisting of six service 
intervals, with two new customers arriving during the first interval and one new customer arriving during 
each of the second, third and fourth intervals. Then if the service discipline is "first-come-first-served" , the 
resulting graph will be 
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whereas if the service discipline is "last-come-first served" , the resulting graph will be 




Nevertheless, the number of vertices and the chromatic number the resulting graph (in this example, 6 and 
3, respectively) are independent of the service discipline. In Section 2, we shall derive asymptotic expansions 
for the moments of these quantities. The leading terms of these expansions are 
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for to > 1, where (2m - 3)!! = (2m - 3) • (2m - 5) • • • 3 • 1 and (-1)!! = 1, 

Ex[if]~lo gT -^, (1.2) 

and, for to > 2, 

where (,{m) = X) n >i V 71 ™ is the Riemann zeta function. It will be noted that Ex[_KT] grows quite slowly 
as A — > 1. If the random variable J denotes the number of customers in the system in equilibrium, then 
Ex [J] = A/(l — A), which grows much more rapidly (see Cohen [C2], p. 181). It may appear paradoxical that 
the maximum number of customers in the system grows more slowly than the mean number of customers, but 
it must be borne in mind that Ex[_K"] is an average over busy periods, whereas Ex [J] is an average over time. 
Indeed, the majority of busy periods have K = 1: after the arrival initiating the busy period, the next event 
determines whether K — 1 (if that event is a service termination) or K > 1 (if that event is another arrival). 
Because A < 1, the former (with probability 1/(1 — A)) is more likely than the latter (with probability 
A/(l - A)). We also note that, since £(2) = tt 2 /6, we have Var^] = Ex[K 2 } - Ex[K} 2 ~ tt 2 /3(1 - A), which 
grows much more rapidly than Ex[_R"] (or even Ex[i^] 2 ). For the chromatic number, corresponding to the 
maximum number of customers in the system simultaneously during the busy period, the results involve 
Lambert series that are generating functions for arithmetical functions arising in number theory: the sums 
of powers of the divisors of positive integers. These are 

Sj(A) = 5>,(n)A n , (1.4) 

n>l 

where cr/(n) denotes the sum of the l-th powers of the divisors of n (see Hardy and Wright [H, p. 239]). 

Our second goal in this paper is to analyze search algorithms connected with stochastic service systems. 
Consider an infinite sequence <Si, S2, ■ ■ ■ of servers in an M /M /oo system. Suppose that each newly arriving 
customer scans this sequence in order and engages the first currently idle server. We are interested in 
the index L of the server Sl engaged by a newly arriving customer in equilibrium. This system has been 
extensively studied by Newell [N2], who suggests that L "is approximately uniformly distributed over the 
interval" [1, A], basing this assertion on the approximation 

Pi[L>l}^ \ l ~ J' if/<A ' (1.5) 

{ 0, if I > A. 

But no error bounds are given for this or other approximations stated by Newell, and not even the fact that 
the first moment has the asymptotic behavior 

Ex[L] ~ ^ (1.6) 



that it would have under the uniform distribution is established rigorously. In Section 3 we shall give a 
rigorous version of (1.5) that will suffice to establish not only (1.6), but also the next term, 

Ex[L] = ^ + ilogA + 0(l), (1.7) 

and more generally 

Ex[L-] = ^ + mAm ~ ll0gA + Q ( A-i) (1.8) 

L J m+ 1 2 v ; v ; 

for to > 1. In particular, we have 

Var[L] = Ex[L 2 ] - Ex[L] 2 

Since the interval [0, 1] is bounded, formula (1.8) shows that the ra-th moment of L/X tends to l/(m + 1) as 
A — > oo for all to > 1, and thus suffices to show that the distribution of L/X tends to the uniform distribution 
on the interval [0, 1]. We note that a problem that is in a sense dual to ours (finding the largest index of 
a busy server, rather than the smallest index of an idle server) has been treated by Coffman, Kadota and 
Shepp [CI]. 

Our final results concern the analogue of the preceding search problem for the M/M/l system. Here 
there is only a single server, but an infinite sequence Wi, W2, ... of waiting stations. A customer arriving 
when the server is busy scans this sequence in order and waits at the first vacant station. When the server 
becomes free and there is at least one customer waiting, it too scans this sequence in order, and serves the 
customer waiting at the first occupied station. We are interested in the index / of the station W/ at which 
a newly arriving customer waits in equilibrium (taken to be zero if the server is idle at the time of arrival). 
(The index of the first station which the server finds occupied (taken to be zero if the service interval initiates 
a busy period) has, of course, the same distribution as I.) We shall show that the distribution of the random 
variable / is closely related to that of the random variable K studied in Section 2, with Ex[/ m ] = XEx[K m ]. 
This fact allows asymptotic expansions for the moments of / to be obtained from those of the moments of 
K in a straightforward way, with the result that the leading terms are the same. 

2. Random Interval Graphs 

Our goal in this section is to determine asymptotic expansions for the moments of the size (number 
of vertices) and chromatic number (number of vertices in the largest clique) for the random interval graph 
corresponding to the M/M/l system. These quantities correspond to the number N of customers served 
during the busy period and the maximum number K of customers in the system simultaneously during the 
busy period. 

The random variable N has a Catalan distribution: 

Pr{N = n} = -^—( 2n - 1 )p n - 1 q n 1 (2.1) 

In — 1 \ n J 

where p = A/(l + A) and q = 1/(1 + A) (see for example Cohen [C2, pp. 190-191], or Riordan [Rl, pp. 64-65]). 
This distribution can be derived as follows. Let J denote the number of customers in the system. When the 
busy period begins, J = 1. During the busy period, J is incremented whenever a new customer arrives, and 



J is decremented whenever a service interval ends and a customer departs, until J = 0, at which time the 
busy period ends (see Figure 1). 




Figure 1. State transition diagram for determining 
the number of customers served during the busy period in M/M/l. 

When J > 1, the probability that the next transition is an arrival is p = A/(l + A), and the probability 
that the next transition is a departure is q = 1/(1 + A). If n customers are served during the busy period, 
there must be n — 1 further arrivals (beyond the one that began the busy period) and n departures, and 
these must occur in such an order that J = for the first time immediately after the last of these n 
departures. The number of such orders is A n — ( ™~ )/(2n — 1) = { n Zi)/n, the n-th Catalan number (see 
for example Comtet [C3, p. 53]). Thus we have (2.1). Since the generating function for the Catalan numbers 
is a(z) = J2 n >i An z n = (l — \/l — 4z)/2, the probability generating function 5at(z) for N is 
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(z) = Y^ Pr[7V = n] . 
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Since for m > 1, we have d m a(z)/dz r 
given by 



_ 1 - VI - Apqz 
2p 

(2m - 3)!!/(l - 4z)( 2m_1 ) /2 , the factorial moments of TV are 



MN(N-l)---(N-m+l)}= 1 - dm W 1 "^-. 



p dz r ' 
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q m 2™- 1 (2m - 3)!! 



VT^Wq (2m - 1] 
_ 2 m - 1 (2m-3)!!A™- 1 
- (1 - A) 2 ™- 1 ' 

because VI -Apq = q-p= (1-A)/(1 + A). Since x m = Eo<;<m { ? } x{x-l) ■ ■ ■ (x-l + 1), where the {7} 
are the Stirling numbers of the second kind, with the generating function X)m>/>o { T } V v. = & ( see 
for example Comtet [C3, pp. 206-207]), and j™} = for m > 1, the ordinary moments of N are given by 



Ex[AH= Yl {7} MN(N-l)---(N-l + l)} 



KKm 



m-i 2'- 1 (2Z-3)!!A'- 1 



E{7}^ 



KKr 



(1-A) 



2Z-1 



Since { m } = 1, this gives the asymptotic formula (1.1) as A — ^ 1. 
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We turn now to the random variable K, the maximum number of customers in the system simultaneously 
during the busy period (counting the customer being served, so that K > 1). The distribution of K is given 

by 

(1 - A) A fc 



PrlK > k] = 



(2.2) 



1 - A fe +! 

(see for example Cohen [C2, pp. 191-193]). This distribution can be derived as follows. Consider a game 
played between two players: P, who begins with v dollars, and Q who begins with w dollars. At each step 
of the game, a biassed coin is tossed; P wins with probability p, in which case Q pays P one dollar, and Q 
wins with the complementary probability q = 1 — p, in which case P pays Q one dollar. The game continues 
until one of the players is ruined (that is, has no money left). It is known that (1) with probability one, 
either P or Q is eventually ruined, and (2), if p ^ q, then the probability that Q is ruined is 



Pr[Q ruined] = 



(q/p) v - 1 
(q/p) v+w - 1 



(2.3) 



(see for example Feller [F, p. 345]). 

Now consider a busy period of the M/M/l queue. The successive events of arrivals and terminations 
of service intervals during the busy period correspond to steps in the game described above. The wealth of 
player P will correspond to the number J of customers in the system, so v = 1. An arrival will correspond 
to a win by player P, so p = A/(l + A), and the termination of a service interval will correspond to a win by 
player Q, so q = 1/(1 + X). Suppose that player Q begins with w = k dollars. Then the event K > k will 
correspond to Q being ruined. Substituting these values in (2.3) yields (2.2) (see Figure 2). 




"~ k ) P > ((K>kj\ 




Figure 2. State transition diagram for determining whether K < k or K > k. 

This correspondence also shows what happens for A > 1. For A = 1 (in which case the busy period is 
finite with probability one, but its expected length is infinite), we have take p = q = 1/2, and have 



This result yields 



so that 



Pr[Q ruined] 



Pr^ >k] = 



v 



v + w 



k + r 

Ex[K] = Y^ p 4 K > k \ 



(2.4) 



/>o 



diverges logarithmically. Of course, for A > 1 (in which case the busy period is infinite with positive 
probability), (2.3) shows that (2.4) diverges linearly. 



Our next goal is to determine the moments of K: 



fe>0 



Writing 

A m (fc) = k m - (k - iy 



. i 

0<Km-l 



E (?)(-i) m - 1 -'* 1 



for the backward differences of the m-th powers, and setting 



rn\n 

T ™( A ) = E i^' ( 2 - 5 ) 



n>l 



summation by parts yields 

Ex[X m ] = E k m Pr[K = k] 

k>0 

= J2 A ™(k + l)Pr[K>k] 

k>0 

A m (k + l)\ k 



fc>0 

A ^ 1 - X> 

1 ~ A V" V" ^ m A I t\m-l-l 3 l X ' 



E £ ?<-i> 



a ^ ^ v i y v y i - a^ 

j>10<Km-l 



= V E (7)(-l) m - 1 -'^(A). (2.6) 

Since Ex[i^ m ] is a linear combination of the 7] (A), it will suffice to determine the asymptotic behavior of 
the sums T m (X). The sums TJ(A) are in fact the Lambert series Si(X) given by (1.4); we have 

-E^E^ 

j>l 4>1 

= EE^ J 
= E^wa", 

n>l 

= Si(X) 

where the inner sum in (2.7) is over integers d dividing n. 



We note that the sums Ti(A) can be expressed in terms of known (albeit exotic) functions of analysis. 
We define the g-gamma function by 

r g (,) = (i- <7 )^n T ^ 

(see for example Gasper and Rahman[Gl, p. 16]). (This function gets its name from the fact that 
lim g _j.i T q (x) — r(x), where T(x) is the Euler gamma function; see for example Whittakcr and Watson 
[W, pp. 235-264].) If wc define the g-digamma function tjjq( x ) as the- logarithmic derivative 

d 

„n+x 
= - log(l - q) + log q Y, 1 _ +x 

of the g-gamma function, then we have 

^ A (l) + log(l-A) 



r (A) 



log A 



To go further, wc define the Z-th g-polygamma function \p q as the l-ih derivative 

^ )(X) =(I^ (X) 
of the g-digamma function. If wc set z = q n+x , then 

d \ fid 

z- 



Since 



dz J \logq dx 






we have 



Summing over n > yields 



J2i~q 



?:>i 



i l q 



log'g \dxj l-q n + 



y^j_ = y i iy q i(n +X ) 

Z-^l \ qIX Z-^i Z-^l * 



1 / d \ l ^ q n+x 



W q \dx) 2-M-grH-* 



1 ^V^(i)+log(l-g) 



log'o v^y io g9 



Thus for / > 1 we have 



>,(0 



T K A)-^« 



log' +1 A' 
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Our next goal is to derive the leading terms (1.2) and (1.3) of the asymptotic expansion of the moments 
of K. To establish (1.2) and (1.3), we begin by deriving 

W^bg^ (2.8) 



Tiw - ;;"\Vi - (2.9) 



and, for I > 1, 

(1-A)' 

Once these formulas are established, it will be clear that the sum in (2.6) is dominated by the term for which 
I = m - 1, so that Ex^" 1 ] ~ m£ m _i(A), and (1.2) and (1.3) follow from (2.8) and (2.9), respectively. 
Our strategy for proving (2.8) and (2.9) will be to approximate the sums Si(X) by integrals 

f°° x l \ x dx 



! 1 - \ x ' 

then then to show that the difference Si(X) — Ii(X) is negligible in comparison with Ii(\). It will be convenient 
to write A = e~ h . The limit A — >• 1 then corresponds to h — > 0. We have 

h = log - 
= log 



l-(l-A) 

1-A. (2.10) 



For I — 0, we have 






/•OO 

/•OO 



/>1 



E 



e 



hi 

i>i 

1 , 1 

- log- 



Substituting (2.10) in this result yields 



h ° 1-A 



W-^log^. (2.11) 



We bound |<So(A) — A)(A)| by the total variation of f(x) = A x /(1 — A 1 ). Since /(a;) decreases monotonically 
from A/(l — A) to as x increases from 1 to oo, we have |So(A) — ^o(A)| < A/(l — A) ~ 1/(1 — A). Since this 
difference is negligible in comparison with (2.11), we obtain (2.8). 



For I > 1, we have 

j x A ctx 



! 1- A x 

{y + l) l \y+ 1 dy 

c v (l\ y* XV+ 1 dy 
^ \i 1-XV+ 1 



r e oyE^ +i) ^ 

E (!)Ef ^ to+1) ^ 



0<i<Z v / j>l ' 

^ \il h l+1 ^ I'+i 

0<i<Z V 7 j>l J 



E ( l ))j^ u ^> ( 2 - 12 ) 



,11 h 

0<i<l 



where Li fc (A) = J2 n>1 X n /n k is the fc-th polylogarithm. Since Lii(A) = log(l/(l — A)) and Li fc (A) — )■ £(fc) as 
A — > 1 for fc > 2, the sum in (2.12) is dominated by the term for which i — I, and we have 



/! Cq + 1) 

*!C(* + 1) 



(i - XY+ 1 



(2.13) 



We bound |£z(A) — /; (A)| by the total variation of f(x) — x l X x /(l — X x ) for < x < oo. As x increases, f(x) 
increases monotonically from to a maximum, then decreases monotonically to 0. Thus the total variation 
of f(x) is twice the maximum. This maximum is 



max f(x) — max 



x l e- h0L 



0<i<oo' 0<:r<co 1 — e 

max 



0<ir<co e' 12 — 1 

1 y* 

— r max 



/V 0<y<oo e y — 1 ' 

Furthermore, y l /(e v - 1) < J!, because e y - 1 = E„>iy"/ n! > V 1 / 11 - Thus l^( A ) - ^( A )l < 
2maxo< I<M /(i) < 2l\/h l ~ 2 Z!/(l — A)'. Since this difference is negligible in comparison with (2.13), 
we obtain (2.9). 

We shall now show how asymptotic expansions, with error terms of the form Of(l — X) R ) for any R, 
can be derived for all of the moments Ex[i4T m ]. The essence of the argument is to use the Euler-Maclaurin 



formula to estimate the difference between Si(X) and Ii(X). This is most conveniently done using a result 
of Zagier [Z]. Indeed, for / > 1, Zagier gives the expansion for Si (A), in terms of the parameter h = — log A 
rather than 1 — A. All that remains for us to do is substitute an expansion for h in terms of 1 — A. For I = 0, 
the expansion for So (A) in terms of h has been given by Egger (nc Endrcs) and Steiner [El, E2], again using 
the result of Zagier. We shall proceed differently, to obtain an expansion involving — log(l — A) rather than 
— log/i. 

Proposition: (Zagier [Z, p. 318]) Let /(x) be analytic at x = 0, with power series f(x) = ^2 r >ob r x r about 
x = 0. Suppose that L°° |/^( x )l dx < oo for all r > 0, where f^ r \x) denotes the r-th derivative of f{x). 
Define F = L f(x) dx. Let g(x) = X) n >i f( nx )- Then g{x) has the asymptotic expansion 

x *-—' (r + 1 

where B r is the r-th Bernoulli number, defined by t/(e t — 1) = ~^2 r>0 B r t r /r\. 
This result is proved by using the Euler-Maclauren formula, 

f f{y)dy= f -f + £ /(n) + ffl+ £ Sf(/ W W-/ (r) (0)) 



l<n<7V-l l<r<i?- 



where -B r (y) is the r-th Bernoulli polynomial, defined by te yt /(e t — 1) = X^r>o ^r(y) t r /r\, and {y} = y — [y\ 
denotes the fractional part of y . (For the Euler-Maclauren formula, the Bernoulli numbers and the Bernoulli 
polynomials, see for example Whittaker and Watson [W, pp. 125-128], where, however, the indexing of the 
numbers and polynomials is different.) The condition L \f^ r '(y)\ dy < oo allows us to let N — > oo, obtaining 



f '/(*)* = £/(»)- £ { ~ ( Z ^r / (r) (°) + (-*)* f / (fl) (i/)^gP^/ 



(r 
If we now write f(xy) instead of /(j/), we obtain 



r.f(xy)dy = j:f(nx)- £ ( ^T / W (0)x f + (-l)V H f {R \xy)^^dy. 



Changing the variable of integration from y to y/x then yields 



irmdy=j2f(nx)- y: h ( 1)r + B iT /(r)(o) " r + (- 1 )***- 1 r f w M B *% /x}) d y- 

The integral on the left-hand side is F, the first sum on the right-hand side is g(x), /^(O) = r! 6 r , and the 
last term on the right-hand side is 0(x R ^ 1 ). Thus 



f = 5 (x)- £ ^h^ + o^i), 

which yields the expansion (2.14). 



x : ' y ' ^ (r + 1) 
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For I > 1, wc define 

m = - 

e- 
Then f(x) is analytic at x = with the Taylor series 



x l 



f> ry.r+1— 1 



r! 

r>0 



and the integral 

X e 



F 



oo „J „- x dx 



1 — <= — x 

i e 



(see for example Whittaker and Watson [W, p. 266]). Furthermore, f( r \x) is a rational function of x and 
e x , in which the degree of the numerator in e x is r, while the denominator is (e x — l) r+1 . Thus f(x) satisfies 
the conditions of the proposition, and we have the asymptotic expansion 

nai + i) , sr (-iy +l - i B r B r+l x r + 1 - 1 

9{x) — + \ WTT) • 

Recalling that A = e~ h , so that h = — log A, we therefore have 

n>l 

J^ (nh) 1 
~ h l ^ e nh - 1 

n>l 

4e/m 

n>l 

ugi + i) , ^ (-ir+'- 1 g r B r+t / t '- 1 

~ /i'+i + ^- r!(r + • ( J 

We note that, if I is odd, then this expansion has only finitely many terms (because B r = for odd r > 3). 
To obtain an asymptotic expansion in terms of f — A, we must substitute the expansion for 1/h: 

1 1 



h — log A 

-i 



log(l - (1 - A)) 
1 -(1-A) 



1 - A l g(l - (f - A)) 



r>0 



where C r is the r-th Bernoulli number of the second kind, defined by t/ log(l + t) = X)r>o ^r t r /r\ (see for 
example Roman [R2, p. 116]). (These numbers are also called the Cauchy numbers of the first kind, and are 
given by C r = L x(x — 1) • • • (x — r + 1) dx; see for example Comtet [C2, pp. 293-294].) 
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For I = 0, we must proceed differently, because 



e x - 1 
has a pole at x — 0. We define 

/*(£) = /(*)_£-! 

1 er 



e x — 1 x 
Then f*(x) is analytic at x — with the Taylor series 



J W ~ ^ (r + l)! 



r>0 

and the integral 

J \e*-l x J 

= 1 

(see for example Whittaker and Watson [W, p. 246]). Furthermore, f* tyR \x) is a rational function of x and 
e x , in which the degree of the numerator in e x is R, while the denominator is Ue x — 1) x) . Thus f*(x) 
satisfies the conditions of the proposition, and we have the asymptotic expansion 

7 ^ {-l) r B r+1 (B r+1 - (-l) r+1 ) x r 
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■w~J+E 



x ^ (r + l) (r+l)! 

r>0 \ < J \ < J 



We therefore have 



W) = E T 



e -nh 



p—nh 
n>l 



n>l 

g— nh | 

^ nh 2-^ e nh - 1 



— nh 1 p —nh 



nh 



n>\ 

~ 1 log -J- + 1 + V ( ~ 1)r * r+1 (i?r+1 ~ ( ' ir+1) *" ■ (2.17) 

/j S l-A/i^ (r + l) (r + l)! l ' 

To obtain asymptotic expansions for the moments of K, we substitute (2.16) into (2.15) and (2.17), 
then substitute the results into (2.6), using the expansion 

1-A 1-A 



A 1 - (1 - A) 

= E( 1 ~ A ) r - 
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Retaining only terms that do not vanish as A — > 1, we obtain 

Ex[X] = log — L- + 7 + O ( (1 - A) log ■ 



1-A ' V 1-A 

and 

MK 2 } = y ^ + log ^ + ( 7 - 1) + O ((1 - A) log ^ 

for the first two moments. Thus we have 
Vax[K] = Ex[iY 2 ] - Ex[K} 2 



log' z r + (1 - 27) log- -Y + 0[ (1 - A) log- 



3 1-A ° 1-A v " °1-A ' V y ° 1-A 



3. Search Algorithms 

In this section we shall analyze search algorithms for M/M/oo and M/M/l systems. We begin by 

studying the distribution of the random variable L, defined as the index of the first idle server S found by 

an arriving customer in the M/M/oo system. Our goal is to prove (1.8), which gives the first two terms in 

the asymptotic expansions of the moments of L. The key to our results is the probability Pr[L > I], which 

is simply the probability that the first I servers Si , . . . , <S; arc all busy. It is well known that this probability 

is given by the Erlang loss formula 

\ l /l\ 



Pr[L > I] 



£o< fe <*A fe A! 
1 
A' 



where 



V- l\ 

(l-k)l\ k 



a= x, n ,:^ fe (3- 1 ) 



0<fc<Z 

(see for example Newell [N, p. 3]). The sum Di can be expressed as an integral, 



oo , _ i 



(see for example Newell [N, p. 7]), and most of Newell's analysis is based on such a representation. But we 
shall work directly with the expression of Di as the sum in (3.1). We shall partition the values of I into 
two ranges. The first, which we shall call the "body" of the distribution, will be < I < Iq = A — s, where 
s = vA. The second, which we shall call the "tail", will be I > lo- 
We begin with the body. We shall establish the estimate 

Pr[L > l} = (1 - //A) + j^ + O (I) + O (^TT^a) (3-2) 

for I < Iq = A — s, where s = VA. We begin by using the principle of inclusion-exclusion to derive bounds 
on the denominator D[. 



13 



We begin with a lower bound. Since 



i(i-i)---(i-k+i)>i k -i j2 ip fe -' 

yo<j</t-i 

= i k -( k V-\ 



we have 



a-E 



1(1 -!)■■■ (I -k + 1) 



X k 

0<k<l 



k , /fc , , r fc _: 



Efl-il 



' \ A7 A ^ 12/ V A 

0<k<l v 7 0<fe<i 



y >l\ k _l + 0((l/\Y) 



For the first sum we have 

^ ( , - 

,X) 1-l/X 

0<k<l v 7 ' 

We note that the logarithm of (l/Xf has a non-negative second derivative for I > 1. Thus (l/X) 1 assumes 
its maximum in the interval < I < Iq for I = 0, I = 1 or / = Iq. Its values there are 0, 1/A and 
(1 — s/A) A ~ s = (1 — l/\/A) A A < e~^ A+1 , respectively. As A — >• oo, the largest of these values is 1/A, so 
we have 0((l/X) 1 ) = 0(1 /X) for < I < l . Thus the first sum is 

v fl\ k = 1 + Q(1/X) 
^ I A/ 1-i/A ' 

0<k<l v 7 ' 

For the second sum we have 

The logarithm of l 2 (l/X) 1 has a non-negative second derivative for ^ > 3, so an argument similar to that used 
for the first sum shows that 0{l 2 {l/X) 1 ) = 0(1/ X) for < I < l . Thus we have 

oh w w " ^ - l ^ 3 



and the lower bound 

1 + Q(1/A) 1 + Q(1/A) 

1 - Z/A X(l-l/X) 
For an upper bound, we have 



Dl > ^^\,'C \Z \!Z - (3-3) 



i(i-i).. ■(i-k + i)<i k ~[ E iM fc_1 + E y'M fc " 2 

iO<i<fe-l / \0<i<j'<*:-l 



• l -l^- l + \&^ 



h2 



(because Eo<i<i<fc-i V = \\Lo<i<k-iJ) - Eo<i<*-i J 2 J / 2 < (Eo<i<fc-iiJ / 2 = (2) / 2 )- Thus 
we have 



X X ^ \2 \X 2A 2 ^ \2 \X 

0<k<l x 7 0<fc<( V 7 V 7 0<fc<Z 
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For the third sum we have 

ECHi) <-E 






K 2 \X ~ ^ \2 \X 

0<k<l v ' v 7 fe>0 

= o- l 



.(1-VA) 5 
and thus the upper bound 

1 + 0(1/A) 1 + 0(1/X) 



A < — T+r^ - — V^tt4 + O 



1-l/X A(l-VA) 3 \X 2 (l-l/X) 5 

Combining this upper bound with the lower bound (3.3) yields 

D _ 1 + 0(1/ X) 1 + 0(1/A) g f 1 



1-Z/A A(l-;/A) 3 VA 2 (1-VA) 5 

To obtain Pr[L > I], we take the reciprocal of Di: 

Pr[£>;] = fl + Q(VA) l + 0(l/A) +0 



1 — l/A A(l-//A) 3 VA 2 (1-VA) 5 

= ( x + w*)) < x - w ( x - xii^uw + ° {-whw 

= (1 + 0(1/A)) (1 - l/X) (i + ^1^ + o (j^jj^ 

B ( 1 + ^» ( (1 ^ /a) + aTI^ + Ka^w))- 

Observing that 0(1/ X) (1 - J/A) = 0(1/ X) and 0(1/A)/A(1 - l/X) = (9(l/A 2 (l - l/X) 3 ), we obtain (3.2). 
We turn now to the tail. We shall establish the estimate 

Pr[L>l] =0(e~ x X l /l\) (3.4) 

for I > X — s, where s = VA. To obtain an upper bound on Pr[L > I], we obtain a lower bound on D;. We 
have 

D - V U 

n n 

- LA - «] ! A'- LA-*J + " ' + [A - 2sJ ! A'- L^-2 S j ' t 3 ' 5 ) 

because I — [X — s\ > I — (A — s) > by assumption and [A — 2s\ > for all sufficiently large A. There 

are [A — 2s\ — [X — 2s\ + 1 > s terms in the sum (3.5). Furthermore, the smallest of these terms is the 

last, because its denominator contains factors of A where the preceding terms contain factors smaller than 

A. Thus we have 

sl\ 
Di > 



LA-2sJ!A'-LA-2 S j 

und, we shall use tl 

for all n > 1 (because the trapezoidal rule underestimates the integral J ± log x dx of the concave function 
logx). This estimate yields 

slle^- 2s i , N 

Di > , . (3.6) 

" e y/[X-2s\ [X - 2s\ L A - 2s J A'- LA-2sj 
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We have 

e |A-2sJ > e A-2s-l^ 

LA - 2sJ LA-2 S j < ( A _ 2 .s)L A - 2;s J 

= A LA-2 sJ(1 _ 2s/A) LA-2,j 

< AL A - 2s J (1 - 2s/X) x - 2s - 1 

< A LA - 2s J e (-2s/A)(A-2s-l) 
<- A LA-2sJ e -2s+4s 2 /A+l 

< A LA ~ 2sJ e~ 2s+5 
and 

a/LA — 2sJ < s. 

Substituting these bounds into (3.6) yields 

l\e x 

Taking the reciprocal of this bound yields (3.4). 

We shall now use (3.2) and (3.4) to prove (1.8). We write 

A m (z) = r - (i - i) m 
= mi™- 1 + o(r~ 2 ) 

for the backward differences of the m-th powers of I. Then partial summation yields 
Ex[L m ] = J2 r P 4 L = l \ 

l>0 

= 5^A m (0Pr[L>J] 
l>o 



J2 m r- 1 Pr[L >l]+0\J2 l ™~ 2 Pt I l > l \ ( 3 - 7 ) 



;>o \ i>o 

This formula shows that we should evaluate sums of the form 



We shall show that 



u n 


l>0 


Pr[L > I}. 


A" 


.+1 


A™ log A 



(3.7) 



Un = , — TT7 , ,, + —^~ + 0(X n ). (3.8) 

(n+l)(n + 2) 2 

Substitution of this formula into (3.7) will then yield (1.8). 

We shall break the range of summation in (3.7) at lo — A — s, where s — a/A, using (3.2) for < I < lo 
and (3.4) for I > lo- Summing the first term in (3.2), we have 

E * n (i-*M) = T E (Ar-r +1 ) 

0<l<l o Q<1<1 

-H(^ +OW) )-(^ +0(5+1) 

1 ^A(A" +1 -(n+l)A" S ) , „^ n ,\ f\ n + 2 -(n + 2)X n + 1 s , „_, 



0(A n ) - ^— ^ + 0(A n+i ) 



A V V "-+1 / V « + 2 

A n+1 



0(A I; 



(n+l)(n + 2) 
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Summing the second term in (3.2), we have 

^ X-l ~ 2-" k 

0<l<l s<k<\ 

s<k<X v 

= A"log- + 0(A n ) 

= ^ + 0(A»), 

where we have used Xa<fc<n V s = log" + 0(1). Summing the third term in (3.2) of course yields 0(X r ' 
Summing the last term in (3.2), we have 

A V l - -A V ( A ~ fc )" 

^ (A-/) 3 ^ k 3 

0<l<l o y ' s<k<\ 

< A " +1 E p 

s<k<\ 

< A " +1 E^ 

fe>s 

A «+! (A + ( ' 



= 0(X n ), 
where we have used ^2 k>n 1/fc 3 = 2/n 2 + 0(l/n 3 ). Combining these estimates, we obtain 

E <" Pr[X > fl = ( „ + 1)(w + 2) + — ^ 4- 0(A»). (3.9) 

Finally, summing (3.4) we have 

s -^l n e~ x X l s ^l n e~ x X l 

l>lo ' l>0 

= 0(X n ), 

because the summation on the right-hand side is the n-th moment of a Poisson random variable with mean 
A, which is a polynomial of degree n in A. Thus 



5^i n Pr[L>l] =0(X n ). 



l>lo 

Combining this estimate with (3.9) yields (3.8) and completes the proof of (1.8). 

Our final goal is to study the distribution of the random variable /, defined as the index of the first 
vacant waiting station W/ found by a customer in the M/M/l system, when the server, upon becoming free 
when at least one customer is waiting, serves the customer at the first occupied waiting station. We shall 
show that the distribution of / is given by 

Pr[/>.] = ^|-^. (3.10) 
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The event I > i is simply the event that the first i waiting stations Wi, . . . ,W{ are all occupied; we have 
1 = when the server is idle. Letting the random variable J denote the number of customers in the system, 
as before, we observe that the event J > i is necessary for the event I > i: if stations Wi, . . . , Wi are 
occupied and the server is busy, there are at least i + 1 customers in the system. Thus we can write 

Pr[J >i]=^2 Pr[/ >i\J = j] Pr[J = j]. 

We shall show that 

Pr[J >i\J = j} = - [_ ~f +1 (3.11) 

for all j > i. Since Pr[J > i] — A* +1 , it will then follow that 

j>i 

1- A 
= T-A^T Pr[J>z] 
_ (1 - A) A' +1 

1 - A*+! ' 

confirming (3.10) 

Before proving (3.11) in the general case, it will be helpful to consider two special cases. If % = 0, 
then the event J > i is sufficient as well as necessary for the event I > i: if the server is busy, an arriving 
customer must wait. Thus Pr[7 > i \ J = j] = 1 for all j > 0, confirming (3.11) in this case. For i = 1, we 
assume that J = j > 1 and ask for the conditional probability that W\ is occupied. If the current arrival 
occurs at time to, we consider the latest transition in the embedded Markov chain for J that precedes to- 
Suppose this previous transition occurs at time t\. If this previous transition was an arrival, then W\ will 
be occupied by it (if it was not already occupied) and thus will be occupied at to- If on the other hand this 
previous transition was a departure, then W\ will be vacated by it (if it was not already vacant) and thus 
will be vacant at to- Thus we must determine the probability that this previous transition was an arrival. 
We claim that the previous transition was an arrival is q = 1/(1 + A) and the probability that it was a 
departure is p = A/(l + A). To prove this claim, we note that the Markov chain for J is reversible; that is, if 
a movie is made of its transitions, the movie run backward is stochastically indistinguishable from the movie 
run forward. (Reversibility follows from the fact that in this Markov chain, transitions occur only between 
adjacent states; that is, J is incremented by an arrival and decremented by a departure.) The previous 
transition was an arrival if and only if it appears as a departure when the movie is run backward, and by 
reversibility this event occurs with probability q = 1/(1 + A) provide that the transition does not involve 
the state J — 0. (When J = 0, the next transition is an arrival with probability one, rather than with 
probability p = A/(l + A).) But since J > 1 at to, we have J > at t\. This proves our claim. We therefore 
have Pr[7 > i \ J = j] = q = 1/(1 + A) = (1 - A)/(l - A 2 ), again confirming (3.11). 

We are now ready to prove (3.11) in the general case. We assume that J = j > i and ask for the 
conditional probability that Wi, . . . , Wi are all occupied. To determine whether this event occurs, we shall 
again trace backward in time through the transitions preceding the current arrival at time to- In this case 
we may have to trace back through arbitrarily many transitions. As we trace backward, we keep track of the 
difference d(t) between the number of arrivals and the number of departures in the interval [t, to). We shall 
continue tracing as long as —1 < d(t) < i, stopping at the latest time t\ such that d(t\) = —1 or d(ti) = i. 
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First, we claim that if d{t\) = i, then Wi, . . . ,Wi are all occupied at time to- To prove this claim, we 
match each departure in [ti,£o) with a later arrival "like parentheses". That is, we associate each departure 
in this interval with a left parenthesis and each arrival with a right parenthesis. Since d(t) > for t\ < t <to, 
there are at least as many right parentheses as left parentheses in any suffix of the resulting string, so all 
the left parentheses can be matched to right parentheses, leaving d(t\) = i right parentheses unmatched. 
We thus have i unmatched arrivals. These arrivals occupy the stations W\, . . . , Wi (if they were not already 
occupied), and any stations that are vacated by subsequent departures are reoccupied by the matching 
arrivals, so these stations all remain occupied at time to- This proves our first claim. 

Next, we claim that if d(ti) = — 1, then at least one of the stations Wi, . . . , Wj is vacant at time to- 
To prove this claim, we define m = max{d(t);ti < t < to}. We have m < i. We define t?, — max{t : t\ < 
t < to and d(t) = m}. Since d{t) < m for t\ < t < t2, we can match arrivals with later departures in this 
interval, leaving at least m + 1 departures unmatched. These departures will vacate stations Wi, . . . , W m +i 
(if they were not already vacant), and any of these stations that are occupied by subsequent arrivals in this 
interval will be revacated by the matching departures, so that these stations will all be vacant at time ti- For 
each of these m + 1 stations Wj, let e(j) denote the difference between the number of arrivals that occupy 
Wj and the number of departures that vacate Wj. We have < e(j) < 1, because arrivals that occupy a 
given station alternate with departures that vacate it. Since dfe) — m, we have Ei<j< m +i e U) = ra - ^ 
follows that e(j) = for at least one value of j, and so Wj remains vacant at time to for this value of j. 
This proves our second claim. 

Finally, we claim that as we trace backward, the probability that the next transition is an arrival is 
always q = 1/(1 + A), and the probability that the next transition is a departure is always p = A/(l + A). 
To prove this claim, we observe that there are at least i + 1 customers in the system at time to, and thus at 
least i + 1 — d(t) customers in the system at any time t for t\ < t < to. Furthermore i + 1 — d(t) > 2, because 
d{t) < i for t\ < t < to- Thus there are at least two customers in the system at every time t e (ti,t ). It 
follows that as we trace backward, none of the transitions we encounter involve the state J = 0. Thus as 
we trace backward, the probability that the next transition is an arrival is always q = 1/(1 + A), and the 
probability that the next transition is a departure is always p = A/(l + A). This proves our third claim. 

These three claims show that the process of determining whether I > i (given that J = j > i) is as 
shown in Figure 3. 




Figure 3. State transition diagram for determining whether I < i or I > i (given that J = j > i). 

Comparison with Figure 2 shows that this process is the same at the process of determining whether 
K > k, except that the roles of p and q are exchanged. This exchange is equivalent to the substitution of 
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1/A for A, so Pr[7 > i \ J = j] is obtained by making this substitution in the expression (2.2) for Pr[K > k] 

(i-i/\)(i/\y 



Pr[7 >i\J = j} = 



(1 - (l/A)*+i) 
1-A 



1 - A«+! ' 

This again confirms (3.11), which completes the proof of (3.10). It follows that the moments of I differ from 
those of if by a factor of A = 1 — (1 — A). This fact allows asymptotic expansions for the moments of I to 
be obtained from those of K, with the result that the leading terms are the same. 

4. Conclusion 

We observe that the search algorithm described in the preceding section for the M/M/l system defines 
a deterministic service discipline distinct from both first-come-first-served and last-come-first-served. It 
would be of interest to determine the distribution of the waiting time W experienced by a customer for this 
discipline, or even the moments of this distribution. We hope to address this question in a future paper. 
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